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Using molecular dynamics simulations, we determine the structure of neutron star crust made 
of rapid proton capture nucleosynthesis material. We find a regular body centered cubic lattice, 
even with the large number of impurities that are present. Low charge Z impurities tend to occupy 
interstitial positions, while high Z impurities tend to occupy substitutional lattice sites. We find 
strong attractive correlations between low Z impurities that could significantly increase the rate of 
pycnonuclear (density driven) nuclear reactions. The thermal conductivity is significantly reduced 
by electron impurity scattering. Our results will be used in future work to study the effects of 
impurities on mechanical properties such as the shear modulus and breaking strain. 
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I. INTRODUCTION 

Neutron stars, collapsed objects half again more mas- 
sive than the sun, are thought to have solid crusts about 
a kilometer thick. Detailed properties of this crust are 
important for many X-ray, radio, and gravitational wave 
observations. Because of the great densities, the elec- 
tronic structure of the crust is likely very simple, con- 
sisting of an extremely degenerate relativist ic Fermi gas. 
The system can be modeled as nearly classical ions inter- 
acting via screened Coulomb, or Yukawa, interactions, 
where the screening length A depends on the electron 
density. 

Indeed many condensed matter systems can be mod- 
eled with Yukawa interactions and much is known about 
the properties of a single component Yukawa system. See 
for example ref. [1 . Because the ion- ion interaction 
is purely repulsive, there is no liquid-gas phase transi- 
tion. However there is a liquid solid phase transition at 
a melting temperature that depends primarily on the 
Coulomb parameter F, 



2^2 



F = 



aT 



(1) 



Here the ions have charge Z, T is the temperature and 
the ion sphere radius a is 



47rn 



1/3 



(2) 



with n the ion density. The system melts at a tempera- 
ture for which F ^ 175 (assuming the screening length is 
relatively large). 
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Neutron stars, that accrete material from a binary 
companion, form new crust from the ashes of nuclear re- 
actions. Simulations of rapid proton capture nucleosyn- 
thesis [2 [3 find a complex composition with many dif- 
ferent ion species. These species then undergo electron 
capture as the material is buried by further accretion to 
greater densities [4 . This still leaves a complex com- 
position of very neutron rich isotopes with many differ- 
ent chemical elements. In this paper, we investigate the 
structure of the resulting solid crust when this complex 
mixture freezes. 

Monte Carlo simulations [5j of the freezing of a classical 
one component plasma (OCP) indicate that it can freeze 
into imperfect body centered cubic (bcc) or face-centered 
cubic (fee) microcrystals. Unfortunetly not much has 
been published on the freezing of a multi-component 
plasma (MCP), although Wunsch et al. study the struc- 
ture of a MCP liquid [6] . There are many possibilities for 
the state of a cold MCP [7^ . It can be a regular MCP lat- 
tice; or microcrystals; or an amorphous, uniformly mixed 
structure; or a lattice of one phase with random admix- 
ture of other ions; or even an ensemble of phase separated 
domains. 

One possibility is that impurities could become frozen 
into random configurations that only relax on very long 
time scales. This could lead to the formation of a glass. 
For example, the binary Fennard- Jones system, with two 
species of different sizes, can form a glass [8^. Here the 
hard core of the interaction keeps the different species 
from diffusing. However the screened coulomb interac- 
tion has a relatively soft 1/r core. This may allow impu- 
rities to diffuse and prevent the formation of a glass. 

In this paper we use molecular dynamics simulations 
to calculate radial distribution functions g{r) and static 
structure factors S{q) to determine the structure of the 
crust. In previous work we determined the chemical sep- 
aration that takes place as the crust freezes. We found 
that the liquid phase is greatly enriched in low charge Z 
ions while the solid is enriched in high Z ions [9 . The 
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distribution of low Z (impurity) ions in the crust may be 
important for the rate of strongly screened thermonuclear 
or pycnonuclear (density driven) reactions. In ref. [10] 
we found that fusion of ^^O + ^^O could be an important 
heat source in the crust. 

In ref. [11 we calculated the thermal conductivity of 
the crust. We found the crust to be a regular crystal 
with a relatively high thermal conductivity. Recently 
the cooling of two neutron stars has been observed after 
extended outbursts \^2^ il3j. These outbursts heat the 
stars' crusts out of equilibrium and then the cooling time 
is measured as the crusts return to equilibrium. The 
surface temperature of the neutron star in KS 1731-260 
decreased with an exponential time scale of 325 ± 100 
days while MXB 1659-29 has a time scale of 505 ± 59 days 
[13 . Comparing these observations, of rapid cooling, to 
calculations by Rutledge et al. [14], Shternin et al. [15], 
and Brown et al. [16 strongly suggest that the crust has 
a high thermal conductivity. 

The shear modulus of a single component system was 
calculated in ref. y/7| where electron screening was found 
to reduce the shear modulus by about 10% compared to a 
pure 1 /r Coulomb system [18] . The shear modulus deter- 
mines the frequency of shear oscillations of the neutron 
star crust. These may have been observed as quasiperi- 
odic oscillations in magnetar giant flares [19]. In the 
future, we will use the results of this paper to calculate 
the effects of impurities on the shear modulus. 

The breaking strain is the deformation of the crust 
when it fails. This determines the maximum height of 
mountains on the surface of neutron stars. These may 
be important sources of gravitational waves for rapidly 
rotating stars [20 . The breaking strain may also be im- 
portant for star quake models of magnetar giant flares 
[21]. In ref. ^22j we examine the effects of impurities on 
the breaking strain. 

The present paper is similar to ref. [11 , however here 
we use a significantly different composition. Ref. pTl 
had so many low Z impurities that phase separation was 
apparently taking place. Therefore, ref. [11 results may 
not be directly applicable to large uniform systems. In 
the present paper we consider a system with fewer impu- 
rities, as described in Section|TI| that may form a uniform 
system. 

Impurities can limit the thermal conductivity. If the 
impurities are weakly correlated then their effect on the 
thermal conductivity can be described by an impurity 
parameter Q [23 , 

Q = (AZ)2= (3) 

This depends on the dispersion in the charge Z of each 
ion. The rp process ash composition of ref. [4 and ref. 
[9j has a relatively large value of Q = 38.9. In this pa- 
per, the composition we use has a smaller, but still sig- 
nificant, value Q = 22.54. This lower value is because 
of a reduction in low Z impurities. Impurity scattering 
can be important at low temperatures where there is less 
scattering from thermal fluctuations. Note that ref. [23] 



assumes the impurities are weakly correlated. If there are 
important correlations among the impurities, for exam- 
ple if there is a tendency for low Z ions to cluster together 
instead of being distributed at random through out the 
lattice, then the effects of impurities on the thermal con- 
ductivity could be different from what is calculated in 
ref. [23 . In this paper we perform MD simulations to 
study the distribution of impurities and their effect on 
the conductivity. 

In section [IT] we describe our molecular dynamics simu- 
lations. Results for the radial distribution function g{r), 
the static structure factor S{q)^ and the thermal conduc- 
tivity are presented in section [lllj We conclude in section 

HYl 



II. MOLECULAR DYNAMICS SIMULATIONS 

In this section we describe our classical molecular dy- 
namics simulations. We begin with a discussion of the 
composition. Schatz et al. have calculated the rapid pro- 
ton capture (rp) process of hydrogen burning on the sur- 
face of an accreting neutron star [2 , see also [3 . This pro- 
duces a variety of nuclei up to mass A ^ 100. Gupta et 
al. then calculate how the composition of this rp process 
ash evolves, because of electron capture and light particle 
reactions, as the material is buried by further accretion. 
Their final composition, at a density of 2.16 x lO^-*^ g/cm^, 
has forty % of the ions with atomic number Z = 34, while 
an additional 10% have Z = 33. The remaining 50% have 
a range of lower Z from 8 to 32. In particular about 3% 
is ^^O and 1% ^^Ne. This Gupta et al. composition is 
listed in the mixture column of Table I in ref. [9 and was 
used for the simulations in ref. [11 . In general, nuclei at 
this depth in the crust are expected to be neutron rich 
because of electron capture. 

Material accretes into a liquid ocean. As the den- 
sity increases near the bottom of the ocean, the material 
freezes. However we found chemical separation when the 
complex rp ash mixture freezes [9|. The ocean is greatly 
enriched in low Z elements compared to the newly formed 
solid. What does chemical separation mean for the struc- 
ture of the crust? Here we make a very simple assump- 
tion and use the composition of the solid phase that was 
found in ref. [9 . [Note that for simplicity we drop chem- 
ical elements with number fraction less than 0.001.] This 
composition is listed in Table [l| and is depleted in low Z 
elements compared to the original Gupta et al. composi- 
tion. For example, we now have only about 1% ^^O com- 
pared to the original 3%. Our composition may not be 
self consistent because we expect chemical separation to 
enrich the ocean in low Z elements and this may change 
the composition of the newly formed crust. This should 
be investigated in future work. 

The electrons form a very degenerate relativistic elec- 
tron gas that slightly screens the interaction between 
ions. We assume the potential Vij{r) between the ith 
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TABLE I: Composition of rapid proton capture nucleosyn- 
thesis ash MD simulations: number fraction Xi of chemical 
element with atomic number Z and mass number A. 
Z A Xi 

8 24 0.0093 

10 28 0.0023 

20 62 0.0023 

22 66 0.0625 

24 74 0.0625 

26 76 0.1019 

27 77 0.0023 

28 80 0.0741 
30 90 0.0949 

32 96 0.0139 

33 99 0.1389 

34 102 0.4306 
36 106 0.0023 
47 109 0.0023 

and jth ion is, 

= ^f^e-'-A% (4) 

where r is the distance between ions and the electron 
screening length is Ag = 7r"'^/^/[2e(37r^ne)"'^/^]. Here ng is 
the electron density. Note that we do not expect our re- 
sults to be very sensitive to the electron screening length. 
For example, the OCP melting point that we found in 
ref. |9j, using a finite Ag, agrees well with the result for 
Ae = oo. 

To characterize our simulations , we define an average 
Coulomb coupling parameter F for the MCP, 

P ^ 

aT ' 

where the ion sphere radius is a = (3/47rn)^/^ and n = 
ne/{Z) is the ion density. The OCP freezes near F = 175. 
In ref. [9 we found that the impurities in our MCP 
lowered the melting temperature until F = 247. Finally, 
we can measure time in our simulation in units of one 
over an average plasma frequency cjp, 

= iL M, ) ' 

j ^ 

where Mj is the average mass of ions with charge Zj and 
abundance Xj (by number). 

We start from initial conditions where we try and min- 
imize arbitrary assumptions about the distribution of im- 
purities (low Z ions) in the solid. A small 432 ion system, 
with composition from Table [l| is started from random 
positions at a high temperature and a relatively high ref- 
erence density n = 7.18 x 10~^ fm~^. Results can be 
scaled to other densities at constant F, Eq. [5] The sys- 
tem is then cooled until it is observed to freeze. Note 



that it is straight forward to crystalize such a small sys- 
tem. Next eight copies of the 432 ion solid are assembled 
into a 3456 ion configuration and this is evolved for a 
short time. Finally eight copies of this 3456 ion system 
are assembled into the final 27648 ion system. The con- 
figuration of this system is shown in Fig. [T] The wavy 
planes of ions in Fig. [l] show that the crystal is strained 
and it is not in equilibrium. 




FIG. 1: (Color on line) Initial Configuration of the 27648 ion 
mixture as described in the text. 



We anneal this 27648 ion system by evolving it accord- 
ing to the cooling schedule in Fig. [2j The simulations 
were performed on a special purpose MDGRAPE-2 board 
[24^ provided by Indiana University's High Performance 
Computing group and took 46 days. The system is first 
heated to near the melting point and then cooled to near 
zero temperature. The final configuration is shown in 
Figs. [3] and [4] We see that the system forms a regular 
crystal with a large distribution of impurities. For exam- 
ple. Fig. |5] shows the configuration of only the Oxygen 
ions. These are seen to be distributed throughout the 
simulation volume. However, there are strong correla- 
tions between the ions that will be discussed in the next 
section. 



III. RESULTS 

We now present results for the radial distribution func- 
tion and static structure factor to characterize the distri- 
bution of impurities in the sample. The radial distribu- 
tion function gij (r) is the probability of finding an ion of 
type j a distance r away from a given ion of type i. It is 
normalized to go to one for large r. We calculate gij by 
histogramming relative distances for 2500 MD configura- 
tions of the 27648 ions where each configuration is sepa- 
rated by a time of 250 fm/c. The dominant species is Se 
{Z = 34), see Table [l| Figure |6] shows the diagonal gii{r) 
for Se-Se correlations. This shows peaks corresponding 
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FIG. 2: Cooling schedule of temperature T versus time t for 
an MD simulation of a 27648 ion system, see text. 




FIG. 3: (Color on line) Final configuration of the 27648 ion 
system at near zero temperature. The medium sized red 
spheres are Oxygen ions, while the small white spheres are 
other below average Z ions, and the large blue spheres are 
above average Z ions. 



to the dominant body centered cubic lattice structure. In 
addition there are dips near r/a = 2.5 and 4 with a the 
mean ion sphere radius. Figure [g] also shows gij{r) for 
correlations between Fe and Se ions. This is very similar 
to that for Se-Se. We conclude that most Fe impuri- 
ties are substitutional and occupy vacant Se lattice sites. 
In contrast, Fig. [6] shows 0-Se and Ne-Se correlation 
functions are significantly different and do not show dips 
near r/a = 2.5 and 4. This suggests that most O and Ne 
impurities are interstitial and occupy positions between 
occupied Se lattice sites. This can be understood if low 
Z impurities are "smaller" than Se ions because reduced 




FIG. 5: Final configuration of only the 256 Oxygen ions (out 
of the total of 27648 ions ) in the system at near zero tem- 
perature. 



Coulomb repulsion allows them to fit into interstitial po- 
sitions. Finally Ca is seen to be an intermediate case. 
Calcium impurities may occupy both substitutional and 
interstitial sites. 

Figure [7| shows diagonal gair) for 0-0, Ti-Ti, Fe-Fe, 
and Zn-Zn correlations as well as Se-Se correlations as 
in Fig. |6] The 0-0 correlation is seen to have a very 
large peak near r/a = 1, note the log scale. This peak 
is consistent with the clustering visible in Fig. [5] This 
suggests that a number of "small" O ions may cluster 
together and take the place of a larger Z ion. To study 
this peak further we show its temperature dependence in 
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FIG. 6: (Color on line) Radial distribution functions gij{r) 
for ions of type i=Se and type j versus r over the mean ion 
sphere radius a at a temperature T = 0.1 MeV. The solid line 
shows g{r) for Se-Se correlations while the other curves show 
correlations between Se and O (dash-dot-dot), Ne (dashed), 
Ca (dot-dashed), and Fe (dotted). 



FIG. 7: (Color on line) Radial distribution functions gii(r) 
for diagonal correlations between ions of type i versus r over 
the mean ion sphere radius a at a temperature T = 0.1 MeV. 
The solid line shows g{r) for Se-Se correlations while the other 
curves show correlations between 0-0 (dash-dot-dot), Ti-Ti 
(dash-dot), Fe-Fe (dashed), and Zn-Zn (dotted). 



Fig. [Sj The peaks in the Se-Se correlation function grow 
sharper with decreasing temperature as the amplitude of 
oscillations is reduced. In contrast the area under the 
first peak in the 0-0 correlation function grows rapidly 
as the temperature decreases. This shows that the O ions 
are becoming very strongly correlated at lower temper- 
atures. This could greatly increase the rate of some py- 
cnonuclear reactions H[l0]. Note that because of these 
strong correlations, it is possible that the O impurities 
may not have fully equilibrated during our simulation. 
Indeed it is even possible that they could phase separate 
at low temperatures, although this could take a very long 
simulation time. 

The static structure factor S{q) describes electron-ion 
scattering. This is important for transport properties 
such as the thermal conductivity, electrical conductivity, 
or shear viscosity. We calculate S{q) directly as a density- 
density correlation function using trajectories from our 
MD simulations, 

5(q) = (p*(q)p(q))-Kp(q))r (7) 
Here the charge density p(q) is, 

1 ^ 7 

with N the number of ions in the simulation and Z^, 
are the charge and location of the ith ion. We evaluate 
the thermal average in Eq. [7| as a time average during 
our MD simulations. We use 2500 configurations that 



are separated in time by 250 fm/c. We also average over 
the direction of the vector q. This calculation is some- 
what time consuming because of the very large number 
of separate q values involved. 

In Fig. |9]we show S{q) at a temperature of 0.1 MeV. 
Although our result has some statistical noise, we see 
peaks that correspond to Bragg scattering from the crys- 
tal lattice. We also show in Fig. [9| as vertical lines of 
arbitrary height, the positions of Bragg peaks expected 
for a body centered cubic lattice of a one component sys- 
tem. The pattern of peaks for our S{q) confirms that 
our lattice is also body centered cubic. However, the 
first peak near qa = 4.2 occurs at a slightly smaller q 
than that for a OOP. This indicates that our unit cell is 
slightly larger and contains slightly more than two ions. 
This is because some of the low Z impurities occupy in- 
terstitial lattice sites. Also shown is a simple fit to one 
component plasma (OOP) results from ref. [25 . Note 
that the OCP fit averages over sharp structures, see for 
example ref. [26]. In addition Fig. [o] shows S{q) for a 
OCP plus the contribution of impurity scattering assum- 
ing the impurities are almost uncorrelated as in ref. [23] . 
This curve assumes an impurity parameter Q = 22.54 
and is below our full S{q) result. This suggests that the 
important correlations that we find between impurities, 
see for example Fig. [8] increase the contribution of im- 
purity scattering to S{q). 

From our S{q) results we calculate the thermal con- 
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FIG. 8: (Color on line) The temperature dependence of the 
radial distribution functions gii{r) versus r over the mean 
ion sphere radius a at temperatures of (bottom to top) T — 
0.3(^reen), 0.2 (red) , and 0.1 (black) MeV. The solid hues 
show g{r) for Se-Se correlations while dashed curves show O- 
O correlations. 



ductivity n as in ref. ^ . 



TTk%TkF 

12a2(Z)A 



(9) 



Here ks is the Boltzmann constant, kp the electron 
Fermi momentum, a the fine structure constant and the 
coulomb logarithm A is. 



A 



2kB 



dq 



ge(g,0)2 



S'{q){l 



(10) 



Here e is the dielectric function of the electrons and S'{q) 
is the inelastic part of S{q)^ see ref. [11 . Our results 
for A and k are listed in Table HIl where we also show 
results (Aimp, /^imp) for the simple fit to OCP results plus 
scattering from uncorrelated impurities. We find that hi 
is somewhat reduced compared to /^imp because of the 
increased contributions of impurity scattering. 



IV. SUMMARY AND CONCLUSIONS 

Using molecular dynamics simulations we have calcu- 
lated the structure of a 27648 ion crystal made from a 
complex composition of rapid proton capture nucleosyn- 
thesis ash that includes many impurities. Even with 
many impurities characterized by a large impurity pa- 
rameter Q = 22.54, we find a regular body centered cu- 
bic crystal with long range order. We do not find an 
amorphous structure. We find that low Z impurities of- 
ten occupy interstitial sites while high Z impurities tend 



FIG. 9: (Color on line) The static structure factor S(q) versus 
momentum transfer q times mean ion sphere radius a. Our 
MD simulation is the black curve. A simple fitting formula 
for a one component plasma is the green dash-dotted line 
while the red dashed line includes impurity scattering from 
uncorrelated impurities. The vertical red dotted lines show 
the positions of Bragg peaks for scattering from a pure one 
component plasma body centered cubic lattice. 



TABLE II: Thermal conductivity k results for a temperature 
of T = 0.043 MeV (5 x lO^K). The original runs were at a ref- 
erence density of 10^^ g/cm^ and the indicated temperatures 
T. The Coulomb logarithm A is defined in Eq. [lO] These 
results have been scaled to the indicated densities. 



A 



Ain 



MeV g/cm^ erg/K cm s erg/K cm s 

0.1 910 0.160 0.121 8.0 x 10^^ 1.57 x 10^^ 2.08 x 10^^ 
0.2 455 0.299 0.253 9.9 x 10^° 4.20 x 10^^ 4.96 x 10^^ 
0.3 304 0.418 0.345 2.9 x 10^° 2.00 x 10^^ 2.42 x 10^^ 



to occupy substitutional sites. There are strong attrac- 
tive short range correlations between low Z impurities 
that grow with decreasing temperature. These correla- 
tions could significantly enhance the rate of pycnonuclear 
(density driven) reactions at high densities. 

The static structure factor S{q) is enhanced by impu- 
rity scattering over S{q) for a one component plasma. 
Furthermore there are important correlations between 
impurities that may invalidate simple models that as- 
sume the impurities are randomly distributed. The ther- 
mal conductivity is reduced by impurity scattering to 
such an extent that if the inner crust is as impure as 
the present simulations for the outer crust show, then 
the electron thermal conductivity in the inner crust will 
be reduced to such an extent that it may disagree with 
interpretations [HI [151 US] of observations [13 of rapid 
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crust cooling. If the interpretation of these observations 
is correct, we conclude that either (a) the large impu- 
rity concentrations in our initial rp ash composition are 
wrong or (b) impurity concentrations are reduced by the 
time material is buried deeper into the inner crust. This 
could be because of nuclear reactions. 

Finally these results will be used in additional work 
to calculate the impact of impurities on the mechanical 
properties of the crust including the shear modulus and 
the breaking strain. The shear modulus is important for 
crust oscillations and the breaking strain may be impor- 



tant for crust breaking models of magnetar giant flares 
and for the stability of mountains that may radiate grav- 
itational waves from rapidly rotating neutron stars. 
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